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We  introduce  a  new  class  of  methods  for  the  Cauchy  problem  for  ordinary  differential 
equations  (ODEs).  We  begin  by  converting  the  original  ODE  into  the  corresponding 
Picard  equation  and  apply  a  deferred  correction  procedure  in  the  integral  formulation, 
driven  by  either  the  explicit  or  the  implicit  Euler  marching  scheme.  The  approach 
results  in  algorithms  of  essentially  axbitraxy  order  accuracy  for  both  non-stiff  and 
stiff  problems;  their  performance  is  illustrated  with  several  numerical  examples.  For 
non-stiff  problems,  the  stability  behavior  of  the  obtained  explicit  schemes  is  very 
satisfactory  and  algorithms  with  orders  between  8  and  20  should  be  competitive  with 
the  best  existing  ones.  In  our  preliminary  experiments  with  stiff  problems,  a  simple 
adaptive  implementation  of  the  method  demonstrates  performance  comparable  to 
that  of  a  state-of-the-art  extrapolation  code  (at  least,  at  moderate  to  high  precision). 
Deferred  correction  approach  based  on  the  Picaxd  equation  appears  to  be  a  promising 
candidate  for  further  investigation. 
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Spectral  Deferred  Correction  Methods  for  Ordinary 

Differential  Equations 


1  Introduction 

The  construction  of  efl&cient,  stable  and  high  order  methods  for  solving  initial  value  problems 
governed  by  systems  of  ordinary  differential  equations  (ODEs)  is,  in  many  respects,  a  mature 
subject  [1,  2,  6,  9,  10,  13].  Existing  methods  for  such  problems  can  be  classified,  coarsely 
speaking,  into  two  groups.  The  first  group  consists  of  intrinsically  high-order  discretization 
schemes  (Runge-Kutta  methods,  linear  multistep  methods,  etc.)  and  the  second  group  con¬ 
sists  of  methods  based  on  accelerating  the  convergence  of  low-order  schemes  through  the 
use  of  Richardson  extrapolation  or  deferred  correction.  For  non-stiff  problems,  there  exist 
extremely  effective  discretizations  of  order  up  to  twelve  or  so,  at  which  point  the  stability 
constraints  imposed  on  the  schemes  become  too  severe.  From  that  point  on,  most  prac¬ 
titioners  recommend  extrapolation.  For  stiff  problems,  the  situation  is  considerably  more 
complicated.  Implicit  Runge-Kutta  methods  possess  excellent  stability  properties,  but  are 
very  expensive  when  high  order  accuracy  is  required.  Implicit  multistep  algorithms  can  have 
very  high  order  convergence,  but  tend  to  have  relatively  poor  stability  properties.  Most 
practitioners,  therefore,  recommend  some  form  of  Runge-Kutta  (or  backward  differentia¬ 
tion)  method  for  orders  up  to  five  or  so,  and  again  turn  to  extrapolation  when  higher  order 
accuracy  is  needed.  These  extrapolation  methods,  while  effective,  are  stiU  expensive,  since 
they  require  computing  a  sequence  of  solutions  on  finer  and  finer  grids.  Although  deferred 
correction  approaches  also  require  computing  a  sequence  of  solutions,  they  are  more  efficient 
in  theory;  the  convergence  rate  can  be  made  to  increase  more  rapidly  and  the  same  underly¬ 
ing  grid  is  used  on  each  sweep.  Because  of  various  numerical  instabilities,  however,  their  use 
has  generally  been  limited  to  the  conversion  of  second-order  accurate  solutions  into  fourth 
or  sixth-order  accurate  ones. 

In  this  paper,  we  present  a  new  version  of  the  deferred  correction  approach.  It  is  based  on 
replacing  the  original  ODE  with  the  corresponding  Picard  integral  equation  and  discretizing 
the  interval  on  which  the  ODE  is  to  be  solved  into  a  composite  Gauss-Legendre  grid.  We 
then  solve  the  integral  equation  approximately  with  either  the  explicit  Euler  method  (for 
non-stiff  problems)  or  the  implicit  Euler  method  (for  stiff  problems)  and  correct  the  solution 
to  higher  and  higher  order  accuracy  by  solving  a  sequence  of  “error”  equations  on  the  same 
grid  with  the  same  marching  scheme.  Because  we  use  spectral  integration  [7,  8],  we  refer  to 
this  class  of  methods  as  spectral  deferred  correction  (SDC)  methods.  For  non-stiff  problems, 
the  approach  results  in  algorithms  of  essentially  arbitrary  order  accuracy.  Moreover,  as 
can  be  seen  in  section  5.1  below,  the  stability  behavior  of  the  resulting  schemes  is  very 
satisfactory.  Our  preliminary  tests  indicate  that  the  schemes  with  orders  between  8  and  20 
are  roughly  competitive  with  the  best  existing  ones.  Our  principal  goal,  however,  is  the  stiff 
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case,  especially  in  the  environment  where  high  precision  is  required.  Since  our  stiff  schemes 
are  driven  by  the  implicit  Euler  method,  we  do  have  to  solve  systems  of  (generally)  non-linear 
equations  at  each  time  step;  unlike  general  implicit  Runge-Kutta  techniques,  however,  we 
do  not  need  to  solve  systems  of  equations  whose  dimensionality  is  greater  than  that  of  the 
underlying  ODE. 

In  certain  respects,  this  paper  should  be  viewed  as  an  experimental  one.  While  the 
convergence  rates  of  the  techniques  we  present  are  easily  proven,  their  stability  properties 
are  established  numerically.  For  orders  up  to  5,  we  have  obtained  schemes  that  are  both 
X-stable  and  A-stable.  For  orders  up  to  30,  we  have  obtained  schemes  that  are  X-stable  and 
A(a:)-stable,  with  a  extremely  close  to  90°  (see  section  5.2  below).  We  have  no  analytical 
reason  to  believe  that  A-stable  schemes  of  order  6  or  higher  do  not,  in  fact,  exist. 

To  fix  notation,  we  assume  that  the  initial  value  problem  to  be  solved  is  in  the  standard 
form 


cp'it)  =  F{tMt))  te[a,b]  (1) 

=  (t>a  ,  (2) 

where  (pa,<p{i)  G  C*"  and  F  :  R  x  C”  ^  Requiring  that  F  G  ^^(R  x  C")  is,  of  course, 
sufficient  to  guarantee  local  existence  and  uniqueness  of  the  problem  (1),  (2).  Since  we 
are  interested  in  high-order  methods,  however,  we  suppose  throughout  that  F  is  sufficiently 
smooth.  Unless  otherwise  stated,  we  assume  that  the  dimension  of  the  system  n  =  1,  since 
it  makes  much  of  the  discussion  less  cumbersome. 

The  structure  of  this  paper  is  as  follows.  In  section  2,  we  introduce  several  analytical  and 
numerical  prerequisites,  in  section  3,  we  describe  the  classical  deferred  correction  scheme 
and  the  difficulties  it  encounters,  and  in  section  4,  we  describe  the  new  spectral  deferred 
correction  approach.  In  section  5,  we  investigate  the  stability  and  accuracy  properties  of 
a  variety  of  SDC  schemes  and  in  section  7,  we  illustrate  the  performance  of  these  schemes 
with  several  examples.  In  section  8,  we  discuss  possible  extensions  and  generalizations  of 
the  approach. 

2  Analytical  and  Numerical  Preliminaries 

In  this  section,  we  summarize  several  well-known  facts  from  numerical  analysis.  First,  sup¬ 
pose  that  we  are  solving  the  problem  (1),  (2)  numerically  on  the  interval  [a,  5],  obtaining  an 
approximate  solution  ^p{h).  The  two  critical  characteristics  of  such  a  scheme  with  which  we 
are  concerned  here  are  its  order  of  accuracy  and  its  (stiff)  stability.  A  numerical  method  is 
said  to  be  of  order  of  accuracy  or  order  k  if,  for  any  sufficiently  smooth  F,  there  exists  a  real 
constant  F  >  0  such  that 


|l^(6)-V^(6)||<F.(6-a)^+U 


(3) 
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The  suitability  of  a  numerical  method  for  stiff  problems  is  generally  analyzed  by  applying  it 
to  the  equation 


ip\t)  =  A  •  <p(t)  t  €  [0, 1] 

9?(0)  =  1  ,  (4) 

The  amplification  factor,  Am{X),  for  A  €  C  is  defined  by  the  formula 

Am{X)  =  ^(1).  (5) 

If,  for  a  given  value  of  A, 

I  Am{X)  |<  1,  (6) 

then  the  numerical  method  is  said  to  be  stable  for  that  value  of  A.  If  a  method  is  stable  for 
all  A  in  the  left-half  plane  (i?e(A)  <  0),  then  the  method  is  said  to  be  A-stable.  A  method 
is  said  to  be  A{a)-stable  if  it  is  stable  for  all  A  such  that  tt  —  a  <  arg{X)  <  tt  -1-  a.  Thus, 
A-stability  is  equivalent  to  A(Q!)-stability  with  ot  =  90®.  Finally,  a  method  is  said  to  be 
L-stable  if 

lim  Am{X)  =■  0.  (7) 

2.1  The  Picard  Integral  Equation 

Integrating  equations  (1)  and  (2)  with  respect  to  t,  we  obtain  the  equivalent  Picard  equation 


ip{t)  =  (paA  [  F{T,(p{T))dT.  (8) 

Ja 

Suppose  now  that  we  have  obtained  an  approximate  solution  ^°(t)  to  (8).  A  measure  of  the 
quality  of  the  approximation  is  given  by  the  residual  function 

S(t)  =  l‘  /(s))&  -  ip\t).  (9) 

Ja 

We  define  the  error  S{t)  by 

6{t)  =  ip{t)  -  (p°{t).  (10) 

Substituting  (10)  back  into  (8),  we  obtain 

/(()  +  S(t)  =  f  F{s,  /(s)  +  S(s))ds,  (11) 

Ja 

or,  after  some  algebraic  manipulation, 

^(^)=  /  [F{s,(p^{s)  +  6{s))- F{s,<p°{s))]ds  +  e{t).  (12) 

Ja 


3 


(13) 


Letting  the  function  G  :  E  x  C  ^  C  be  given  by 

G(t,S)  =  F(t,v\t)  +  «(())  -  F((,  v°(i)), 
we  can  rewrite  (12)  in  the  form 

S{t)-  f  G{s,S{s))ds  =  e{t),  (14) 

J  a 

which  is  a  Picard- type  integral  equation  like  (8). 

2.2  Euler  Methods  for  the  Piccird  Equation 

Suppose  that  to, fi,  ^2,  •  •  • , <m+i  is  a  refinement  of  the  interval  [a,  b]  with 

^0  =  o,  (15) 

tm+l  —  b,  (16) 

to<ti<t2  •••  <tm<tm+l.  (17) 

Then,  the  explicit  Euler  (or  forward  Euler)  method  for  the  solution  of  the  ODE  (1)  or  the 
integral  equation  (8)  is  given  by  the  formula 

9i+i  =  V’i  +  hi  ■  F(ti,  ipi)  (18) 

hi  —  ti,  (19) 

for  i  =  0, 1,  •  •  • ,  m.  The  implicit  (or  backward)  Euler  scheme  for  the  solution  of  (1)  is  given 
by  the  formula 

¥’i+i  =  ^^1  +  ^1'  P’(i,-+i,  V’i+i)-  (20) 

Similarly,  the  explicit  Euler  method  for  the  solution  of  (14)  is  given  by 

=  hi  +  hi  •  G{ti,  Si)  +  (£(fi+i)  -  £(<i))5  (21) 

and  the  implicit  Euler  scheme  for  the  solution  of  (14)  by 

^j+i  =  Si  +  hi  •  G{ti+i,Si+i)  -f  (e(t,-+i)  —  e{ti)).  (22) 

Definition  2.1  Given  the  function  G  :  E  x  C  ^  C  defined  in  (13)  and  the  vector  of  data 
e  =  {e(ti),£(t2),  •  •  •  e(tm)}  €  we  define  the  map  C^x-p  :  G^(E  x  C)  x  — »•  by 

GexpiG^s)  =  S, 

where  S  =  (5i,  ^2,  •  •  •  ?  ^m)  is  the  vector  of  corrections  produced  by  the  scheme  (21 ).  Similarly, 
we  define  the  map  Gimp  :  G^(E  x  C)  x  €”*  — »■  O’"  by 

Gimp{G,e)  =  S, 

where  S  =  {Si,  S2,---,  S^)  is  the  vector  of  corrections  produced  by  the  scheme  (22). 
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A  full  description  of  the  marching  schemes  corresponding  to  Cexp  and  Cimp  requires  that 
we  specify  how  the  values  e(tj)  are  actually  computed.  For  this,  we  will  require  stable  and 
high-order  accurate  methods  for  interpolation  and  integration. 


2.3  Spectral  Integration,  Differentiation,  and  Interpolation 


Given  a  natural  number  m,  we  will  denote  by  ri,r2,-“  ,rm  the  m  Gauss-Legendre  nodes 
on  the  interval  [—1, 1]  (see,  for  example,  [17]).  For  an  interval  [o,  b]  €  R,  we  will  denote  by 
si, 52?  •  •  •  5 the  m  Gaussian  nodes  on  the  interval  [a,  6],  given  by  the  formula 


Si  = 


b  —  a 


Ti  + 


6  d*  £Z 
2  ■ 


(23) 


Suppose  now  that  is  a  strictly  increasing  sequence  of  points  in  R,  and  that 

with  each  point  U  is  associated  a  function  value  (fi.  Let  (^  =  ((^i,  </?2,  • .  •  >  Then,  for  any 
point  t  G  R,  we  will  denote  by  L™  :  C”  x  R  (D  the  usual  Lagrange  interpolant  defined  by 
the  formula 

m 

L‘^{(p,t)  =  Y^Ci{t)^(pi,  (24) 

t=i 

where  the  functions  c,(t)  are  given  by 


ci{t) = n 

j^i 


t-tj 
ti  tj 


(25) 


Definition  2.2  Let  jP  :  R  — >  (D  and  let  the  vector  f  =  {/i,/2,  —  fm}  be  defined  by  the 
formula 

fi  =  F{ti).  (26) 

ly  e  =  {ei,  62,  •  •  • ,  Cto}  is  defined  by 

^  T  {f,ti),  (27) 

then  the  linear  mapping  D™‘  :  (ST'  C*"  for  which 

e  =  D^(f)  (28) 

will  he  referred  to  as  the  differentiation  matrix.  If  g  =  {go,9i,  "•  ,9m}  is  defined  by 

9i  =  £'^L^{f,t)  dt,  (29) 

then  the  linear  mapping  :  C’"  — >  C”  for  which 

S  =  S”)/)  (30) 

will  be  referred  to  as  the  integration  matrix. 
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If  the  function  F  is  a  polynomial  of  degree  m  —  1,  and  the  vector  /  is  defined  as  in  (26), 
then 

Fit)  =  L^(f,t),  (31) 

and  the  operators  and  S”^  are  exact. 

Remark  2.1  The  formulae  (28)  and  (30)  are  not  numerically  stable,  unless  the  points 
are  chosen  with  some  care;  using  equispaced  nodes,  for  example,  leads  to  the 
well-known  Runge  phenomenon.  On  the  other  hand,  with  a  suitable  choice  of  nodes,  the 
operators  defined  in  (28)  and  (30)  become  extremely  effective  numerical  tools.  The  most 
popular  choices  are  the  Chebychev  and  Gauss-Legendre  nodes,  for  which  the  matrices  D'^, 
are  usually  referred  to  as  spectral  differentiation  and  spectral  integration  matiices,  respec¬ 
tively.  We  refer  the  reader  to  [7,  8, 17]  for  a  detailed  discussion  of  their  numerical  properties. 
Here,  we  simply  observe  that  spectral  integration  is  an  extremely  useful  tool;  the  maxirmim 
eigenvalue  of  is  bounded,  and  its  minimum  eigenvalue  is  of  the  order  0(l/m^).  Spec¬ 
tral  differentiation,  while  widely  used,  is  somewhat  limited  by  the  fact  that  the  maximum 
eigenvalue  of  is  of  the  order  0{rv?),  which  renders  the  operator  ill-conditioned. 

Remark  2.2  It  is  easy  to  see  that  the  matrices  are  dense  for  any  distribution  of 

points  ti,  t2,---  tm-  Since  applying  a  dense  mxm  matrix  to  a  vector  requires  operations, 
the  procedure  can  become  quite  expensive.  When  <1,^2, •••im  are  Chebychev  nodes  on 
the  interval  [a,  6],  however,  the  Fast  Fourier  Transform  (FFT)  can  be  used  to  apply  the 
matrices  S”",  D’”  to  arbitrary  vectors  using  0(m  log  m)  operations.  For  more  general  point 
distributions,  a  somewhat  less  efficient  0[m  log  m)  scheme  for  the  application  of  the  matrices 
5™  and  to  arbitrary  vectors  can  be  found  in  [4].  In  the  present  paper,  we  will  be  using 
relatively  short  sequences  (m  «  16),  so  that  the  cost  of  applying  5”"  and  to  arbitrary 
vectors  will  not  be  a  major  issue. 

In  section  6,  we  will  need  one  more  numerical  tool.  Let  ^i,  r2,  •  •  • ,  denote  the  Gauss- 
Legendre  nodes  on  the  interval  [-1, 1].  We  then  define  the  m  x  m  matrix  by  the  formula 

=  (32) 

where  P,-  denotes  the  Legendre  polynomial  of  order  j.  Note  that  maps  the  vector 
(«!,...,  arn)  into  the  vector  (/i,  /2, . . . ,  /m),  where 

m 

fi  =  ^ayP,_i(r,-). 

The  matrix  is  non-singular  [7]  and  we  will  denote  its  inverse  by 

W*"  =  (y-)-L  (33) 

Given  a  polynomial  Q  of  degree  m  -  1,  it  is  clear  that  the  matrix  W'"  maps  the  values  of  Q 
at  the  nodes  ri,  r2,  •  •  • ,  into  the  coefficients  of  its  Legendre  expansion. 


6 


3  Classical  Deferred  Correction 


Suppose  now  that  we  define  a  grid  on  the  interval  [a,  b]  with  (m  +  1)  equispaced  nodes  U 
given  by 

U  =  a  +  i  ‘  h  i  —  0, m  ,  (34) 

where  h  =  [b  —  a)/m  is  the  step  size  and  that  we  wish  to  solve  the  ordinary  differential 
equation  (1),  (2)  on  this  grid.  A  kth.  order  accurate  method  will  yield  an  approximate 
solution  T]  =  with 

Tii  =  ^{ti)  +  0{h’‘)  .  (35) 

This  defines  the  unique  mth  order  polynomial  L”^{r},t)  which  interpolates  the  discrete  ap¬ 
proximate  solution  values  r]i  at  the  designated  grid  points  t,-.  We  can  then  define  an  error 
function 

=  (36) 

which  clearly  satisfies  the  differential  equation 

=  fit.  Sit)  +  X»(,,  <))  -  i)  (37) 

«(0)  =  0.  (38) 

We  can  now  solve  this  equation  for  the  error  function  by  the  same  kth  order  method  as  used 
for  the  original  problem.  In  other  words,  we  generate  a  sequence  of  values 

TTi  fa  S{ti)  i  =  1, ...,  m  ,  (39) 

on  the  same  grid  as  used  previously.  It  is  well-known  that  the  “corrected”  approximation 

rji  TTi  fa  y{ti)  i  =  1,  •  •  • ,  m  (40) 

is  of  {2ky^  order  accuracy  [3,  5,  14,  15,  16,  19]. 

Iterated  deferred  correction  proceeds  by  computing  a  new  polynomial  interpolant  to  the 
updated  approximate  grid  values  (t,-,?7i  -I-  tt,),  defining  a  new  error  function,  and  solving  a 
new  correction  equation  of  the  same  form  as  (37)  above. 


Algorithm:  Deferred  Correction 

Comment  [Compute  initial  approximation] 

Using  a  order  method,  compute  an  approximate  solution  « 
at  the  grid  points  t,-,  i  =  1, ...,  m  on  the  interval  [0,T]. 
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Comment  [Compute  successive  corrections.] 
do  j  =  1,. . .  ,J 

1)  Compute  the  interpolating  polynomial  t). 

2)  Define  the  error  function  8{t)  =  <f{t)  — 

3)  Form  the  error  equation  S'{t)  = 

with  6(0)  =  0. 

Comment  [Note  that  the  values  of  the  derivative  at  the  grid  points  are  contained 

in  the  vector  where  is  the  differentiation  matrix.] 

4)  Using  a  order  method,  compute  an  approximate  solution  tt,  »  6(t,) 

at  the  grid  points  tj  on  the  interval  [0,r]. 

5)  Define  a  new  approximate  solution  +  tt,-  . 

enddo 

At  the  end  of  this  procedure,  the  error  is  of  the  order 

Of  course,  this  process  can  only  be  repeated  so  long  as  and  t)  are 

sufficiently  accurate.  The  usual  estimate  of  the  order  of  accuracy  obtained  with  iterated 
deferred  correction  is  [3] 

0(h™n[(*^+i)-<'.»n])  .  (42) 

There  are  two  independent  factors  which  prevent  the  use  of  large  m,  and  which  have 
prevented  large  numbers  of  iterations  from  being  used  in  practice.  The  first  problem  relates 
to  the  instability  of  approximation  at  equispaced  nodes;  as  mentioned  earlier,  the  procedure 
is  numerically  ill-conditioned  (the  Runge  phenomenon).  The  second  problem  is  that  the 
procedure  involves  numerical  differentiation  in  the  construction  of  the  new  right-hand  side  for 
each  error  equation.  Differentiation  introduces  subtle  instabilities  which  prevent  the  effective 
use  of  large  m  (for  related  phenomena,  see  [11,  18]).  The  difficulty  of  interpolation  is  easily 
eliminated  through  the  use  of  Legendre  polynomials.  The  need  for  numerical  differentiation 
is  eliminated  by  using  the  Picard  equation  as  our  starting  point. 
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4  Spectral  Deferred  Correction  schemes 

Suppose  now  that  we  are  given  an  approximate  solution  9?^®^  on  the  interval  [a,b].  We  have 
already  described  the  error  equation  (14)  which  arises  from  the  Picard  formulation  of  the 
original  ordinary  differential  equation,  and  it  remains  only  to  complete  our  description  of 
the  discretization  process. 

In  the  remainder  of  this  paper,  we  will  use  the  grid  corresponding  to  the 

standard  Gauss-Legendre  nodes  on  [a,b].  will  be  used  to  denote  the  jth  approximate 
solution 

9?b']  =  «  (9?(5i), 97(62),..., v?(s,„)), 

^  will  be  used  to  denote  the  m-vector  (9?a,¥’o? •  •  •  iV’o)?  and  F{ip^^)  will  be  used  to  denote 
the  vector 

(F(si,  9?^'H5i)),  F(s2,  <^^^(52)),  • . .  F{s^,  9s[^'J(s„j))). 

The  residual  function  £(t)  defined  in  (9)  will  be  approximated  by  the  vector  defined 

by  _ 

F(9.t^1)  -  +  Tp-a-  (43) 

Observe  that  (43)  is  obtained  from  (9)  by  replacing  exact  integration  with  spectral  inte¬ 
gration.  We  may  now  proceed  with  the  construct  of  high  order  schemes  for  both  stiff  and 
nonstiff  ODEs. 


Algorithm:  Spectral  Deferred  Correction 
Comment  [Compute  initial  approximation] 

For  nonstiff/stiff  problems,  use  the  forward/backward  Euler  method  to  compute  an  approximate 
solution  (pf^  «  9>(si)  at  the  grid  points  si, . . . ,  on  the  interval  [o,  6]. 

Comment  [Compute  successive  corrections.] 

do  j  1,. « « ,J 

1)  Compute  the  approximate  residual  function  cr(9)b-i]). 

2a)  For  nonstiff  problems,  compute  ^ =  Cexp{G,  cr{<p^~^^))  as  in  section  2.2. 

2b)  For  stiff  problems,  compute  =  Cimp{G,(T{(p^~^^))  as  in  section  2.2. 

3)  Update  the  approximate  solution  tp^  =  4.  ^bl. 

enddo 
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Definition  4.1  For  nonstiff  problems,  the  numerical  method  outlined  in  the  preceding  algo¬ 
rithm  using  m  nodes  and  J  correction  steps  will  be  denoted  by  EuExp^.  The  approximate 
solution  generated  by  the  scheme  will  be  denoted  by  EuExp^{F,(fa).  For  stiff  problems, 
the  numerical  method  outlined  in  the  preceding  algorithm  using  m  nodes  and  J  correction 
steps  will  be  denoted  by  Eulmp^.  The  approximate  solution  gP]  generated  by  the  scheme 
will  be  denoted  by  EuImp^(F,ipa). 

As  in  the  classical  deferred  correction  case,  it  is  straightforward  to  obtain  the  following 
result  [3]. 

Theorem  4.1  For  any  sufficiently  smooth  function  F  :  R  x  C  C  and  any  natural  numbers 
m,k,  each  of  the  approximations  EuExpl^[F,ip a)  and  EuImp'^{F,ipa)  converge  to  the  exact 
solution  with  order  of  accuracy  min(m,  J  +  1). 

Remark  4.1  An  apparent  drawback  of  the  schemes  EuExp^  and  Eulmpi^  as  tools  for  the 
solution  of  the  initial  value  problem  (1),  (2)  on  the  interval  [a,  b]  is  the  fact  that  the  nodes 
*^1 9  •  •  •  5  lie  inside  the  interval,  so  that  no  solution  is  generated  at  the  endpoint  b.  This 
problem  is  easily  remedied  by  using  the  interpolating  polynomial  to  obtain 

<Pb  =  L^{EuExpi{F,  <pa),  b).  (44) 

If  the  solution  is  desired  at  an  arbitrary  point  t  in  the  interval  [a,  6]  we  again  use  the 
Lagrange  interpolant  L'^{EuExp'^{F,ipa.),t). 

Remark  4.2  (Systems  of  DDEs)  When  considering  systems  of  DDEs,  one  simply  per¬ 
forms  the  interpolation  and  integration  operations  componentwise.  The  function  evaluation 
and/or  inversion  which  is  required  at  each  step  of  (21)  or  (22)  is  obviously  more  complicated, 
but  has  no  effect  on  the  overall  structure  of  the  schemes  EuExpf^  and  Eulmpf^. 

Remark  4.3  In  the  stiff  case,  a  general  implementation  of  (22)  will  involve  the  solution  of 
a  nonlinear  equation  (or,  more  generally,  a  system  of  nonlinear  equations).  Normally,  this  is 
done  using  some  form  of  Newton’s  method  (see,  for  example,  [12]).  One  is  then  confronted 
with  numerous  issues  such  as  the  choice  of  the  initial  approximation,  error  control,  iteration 
count,  and  the  frequency  with  which  the  Jacobian  of  G  is  recomputed  and  inverted.  In 
the  context  of  (22),  most  of  these  issues  are  simple  (except  when  evaluating  the  initial 
approximation  Indeed,  since  the  correction  8  is  expected  to  be  small,  0  can  be  chosen 
as  the  initial  approximation,  only  one  iteration  of  the  Newton  procedure  is  required  at  each 
step,  and  the  recomputation  and  inversion  of  the  Jacobian  can  be  bypassed  once  the  accuracy 
of  the  approximation  is  of  the  order  y/e,  where  e  is  the  desired  accuracy  of  the  computation. 
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4.1  Behavior  at  Large  |  A  |  and  Extrapolation 

While  the  methods  EuExp^  and  Eulmp'^  tend  to  be  quite  satisfactory  for  non-stiff  and 
mildly  stiff  problems,  for  strongly  stiff  problems  we  will  need  an  additional  modification.  We 
start  with  the  following  obvious  theorem. 

Theorem  4.2  For  any  pair  of  natural  numbers  m,  J,  the  amplification  factor  Am{X)  asso¬ 
ciated  with  the  scheme  Eulmp^  is  a  rational  function  of  X.  Furthermore,  there  exists  a  real 
number  fi{m,  J)  such  that 

lim  Am{X)  =  p,{m,J).  (45) 

|A|“>-oo 

For  all  combinations  m,  J  we  have  tested,  /f(m,  J)  <  1,  making  them  acceptable  for  stiff 
problems.  We  have  not  encountered  any  combinations  m,  J  for  which  the  scheme  Eulmp^ 
is  //-stable,  though  some  are  A(a)-stable  with  fairly  large  a.  Fortunately,  the  above  theorem 
provides  a  mechanism  for  combining  two  different  schemes  with  different  m,  J  to  obtain  in¬ 
stable  schemes;  for  reasons  the  authors  do  not  completely  understand,  the  resulting  schemes 
also  tend  to  have  much  improved  A(q:)  -  stability. 

Corollary  4.3  Suppose  that  mi,  ji,  m2,  ^2  are  four  positive  integer  numbers  such  that  /u(mi,  ji)  ^ 
M(m2,ji2),  and  EuC is  the  scheme  for  the  solution  of  the  problem  (1),  (2)  defined 
by  the  formula 


Then  EuComlP^'^J^^  is  L-stable. 

Remark  4.4  We  have  not  carried  out  a  systematic  investigation  of  the  properties  of  spec¬ 
tral  deferred  correction  schemes  based  on  Gauss-Radau  or  Gauss-Lobatto  discretization, 
which  includes  one  or  both  endpoints.  We  have,  however,  experimented  with  Chebychev 
discretization  in  this  context,  and  obtained  results  very  similar  to  those  reported  in  this 
paper.  The  highest  order  EuComb  scheme  which  we  found  to  be  A-stable  scheme  appears 
to  be  3,  whereas  with  Gaussian  nodes,  EuCombl’^  is  A-stable  and  has  order  5.  We  have  no 
analytical  results  explaining  this  difference,  and  for  most  practical  purposes,  Gaussian-based 
and  Chebychev-based  schemes  are  very  similar.  Because  of  this  difference,  we  conjecture 
that  there  may  exist  nodes  which  lead  to  A-stable  schemes  of  order  higher  than  5. 


4.2  Composite  Schemes  and  Stability  issues 

When  considering  a  general  purpose  solver  for  the  initial  value  problem  (1),  (2)  on  the 
interval  [a,  6],  it  is  rarely  reasonable  to  use  a  single  global  mesh.  Thus,  we  assume  that  the 
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interval  [a,  6]  is  subdivided  into  a  collection  of  subintervals  [a,-,  6,],  such  that  a,+i  =  bi,  and 
apply  one  of  the  schemes  EuExp,  Eulmp  or  EuComb  (with  a  reasonably  small  m)  on  each. 
Such  an  algorithm  is  very  similar  to  a  Runge-Kutta  method;  it  is  essentially  a  single-step 
algorithm  with  limited  storage  requirements,  it  is  easy  to  implement  adaptively,  and  its 
stability  properties  are  not  obivous  a  priori.  It  turns  out  that  for  non-stiff  problems,  EuExp 
works  extremely  well  for  a  variety  of  combinations  of  the  parameters  m,  J.  For  stiff  problems, 
schemes  of  the  EuExp  type  are  obviously  useless,  since  they  are  driven  by  the  explicit  Euler 
method.  Schemes  based  on  Eulmp  result  in  acceptable  methods  for  certain  choices  of  m,  J. 
Unfortunately,  for  larger  values  of  m,  J,  the  stability  properties  of  Eulmp^  deteriorate 
rapidly.  Finally,  schemes  based  on  EuComb  result  in  acceptable  stability  properties  for 
many  values  of  mi,ji,  m2, 72-  All  such  methods  are  i-stable,  and  many  combinations  of 
result  in  nearly  A-stable  schemes  (see  section  5.3  below).  While  no  general 
analysis  of  such  methods  has  been  carried  out,  our  experiments  appear  to  indicate  that 
there  exist  schemes  of  this  type  that  axe  of  arbitrarily  high  order  and  are  A(a)-stable  with 
a  extremely  close  to  90®. 

5  Stability  and  Accuracy  Properties  of  Selected  Schemes 

We  have  implemented  the  schemes  EuExp,  Eulmp,  and  EuComb  in  FORTRAN,  and  we 
have  conducted  a  number  of  numerical  experiments  with  the  resulting  codes  in  order  to 
elucidate  their  performance.  The  following  terminology  is  used  in  this  section.  The  stability 
region  associated  with  a  numerical  scheme  for  the  solution  of  the  equation  (4)  is  defined  to 
be  the  subset  of  the  complex  plane  C  consisting  of  all  A  such  that  on  the  interval  [0, 1],  the 
amplification  factor  define  in  (5)  satisfies  A.m(A)  <  1.  For  a  given  e  >  0,  the  accuracy  region 
associated  with  a  numerical  scheme  is  defined  to  be  the  subset  of  C  consisting  of  all  A  such 
that,  when  the  scheme  is  applied  to  the  equation  (4)  on  the  interval  [0, 1], 

\  (p{b)  -  <f{b)  \<  e.  (47) 

Since  both  (p{b)  and  y?(&)  are  analytic  functions  of  A,  it  follows  from  the  maximum  principle 
that  both  the  stability  and  accuracy  regions  have  well-defined  boundaries. 

5.1  Stability  and  Accuracy  Properties  of  EuExp  schemes 

We  first  compute  the  boundaries  of  the  stability  and  accuracy  regions  for  the  schemes 
EuExp^  for  several  choices  of  the  parameters  m  and  J  (Figures  1  —  4).  It  should  be  noted 
that  the  stability  regions  are  compact;  in  other  words,  they  are  stable  inside  the  boundaries 
marked  j4m(A)  =  1.  Several  observations  can  be  made  from  these  figures  and  from  the  more 
detailed  numerical  experiments  we  have  performed. 
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1.  In  all  cases,  the  stability  condition  is  dominated  by  the  accuracy  condition.  In  other 
words,  whenever  Re{X)  <  0  and  the  scheme  achieves  a  reasonable  accuracy,  the  scheme  is 
stable. 

2.  The  sizes  of  both  the  stability  and  accuracy  regions  grow  with  the  order  of  the  scheme; 
when  A  is  purely  imaginary,  the  scheme  of  order  20  requires  about  20  nodes  per  wavelength 
to  achieve  11  -  digit  precision. 

Remark  5.1  The  Euler  method  is  obviously  not  the  most  efficient  solver  to  which  the 
deferred  correction  approach  can  be  applied.  We  have  experimented,  for  example,  with 
explicit  Adams  methods  of  orders  up  to  6,  and  have  obtained  improvements  of  up  to  a  factor 
of  three  in  terms  of  the  number  of  function  evaluations  required. 


5.2  Stability  and  Accuracy  Properties  of  Eulmp  schemes 

For  a  number  of  combinations  m,J,-we  have  numerically  constructed  the  boundaries  of  the 
stability  and  accuracy  regions  for  the  schemes  Eulmp'^  (Figures  5  -  12).  The  regions  of 
stability  of  these  schemes  extend  to  infinity;  in  other  words,  they  are  stable  outside  the 
boundaries  marked  Am{X)  =  1.  It  is  worth  noting  that,  in  most  cases,  the  regions  of  insta¬ 
bility  are  very  much  larger  than  the  regions  of  accuracy.  Thus,  for  each  scheme  we  present 
two  figures.  The  first  is  on  a  relatively  coarse  scale,  depicting  the  stability  region.  The 
second  is  on  a  much  finer  scale,  depicting  the  accuracy  regions  for  two  selected  accuracies; 
in  the  latter  case,  the  boundary  of  the  stability  region  is  virtually  indistinguishable  from  the 
imaginary  axis.  Each  of  the  figures  carries  a  legend,  specifying  detailed  stability  characteris¬ 
tics  of  the  scheme  (all  of  the  schemes  Eulmp^  are  j4(Q;)-stable,  and  the  legends  specify  the 
approximate  values  of  a,  obtained  numerically).  None  of  the  schemes  Eulmp^  are  T-stable; 
the  legends  specify  the  value  of  p  for  each  of  the  schemes  (see  (45)). 

Several  observations  can  be  made  from  the  figures  5  -  12. 

1.  There  exist  Eulmp'^  schemes  that  are  A-stable  of  order  up  to  four  (such  as  Eulmp^'). 
The  scheme  Eulmpl  is  A(a)-stable  with  a  >  89.5°;  for  most  practical  purposes,  such 
a  scheme  can  be  viewed  as  A-stable.  For  higher  orders,  the  A  -  stability  properties  of 
Eulmp'^  deteriorate  rapidly  (see  Figures  9-12). 

2.  None  of  the  schemes  Eulmp^  are  T-stable.  However,  in  all  cases  we  have  studied,  p 
is  less  than  1/2;  while  L  -  stability  (corresponding  to  /i  =  0)  is  very  desirable,  p  <lf2 
guarantees  a  rate  of  decay  that  is  sufficient  in  many  cases. 

3.  The  accuracy  regions  for  all  methods  we  have  tested  are  very  satisfactory,  for  both  real 
and  complex  A.  It  is  easy  to  see,  for  example,  from  Figure  8  that  Eulmpl  scheme 
of  order  6)  requires  about  18  nodes  per  wavelength  to  obtain  3  digits;  the  number 
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increases  to  about  40  nodes  per  wavelength  to  obtain  5  digits,  indicating  the  need  for 
a  higher  order  scheme.  Such  schemes  are  discussed  in  the  following  subsection. 

5.3  Stability  and  Accuracy  Properties  for  EuComh  schemes 

For  a  number  of  combinations  mi,m2,  we  constructed  numerically  the  boundaries  of 
the  stability  and  accuracy  regions  for  the  schemes  (Figures  13  —  24).  As  in  the 

case  of  the  EuComb  schemes,  the  stability  regions  of  these  schemes  extend  to  infinity;  they 
are  stable  outside  the  boundaries  marked  Am{X)  =  1.  As  for  the  simpler  Eulmp  schemes, 
the  regions  of  instability  are  generally  very  much  larger  than  the  regions  of  accuracy.  Thus, 
we  again  present  two  figures  for  each  case.  The  first  is  on  a  relatively  coarse  scale,  depicting 
the  stability  region.  The  second  is  on  a  much  finer  scale,  depicting  the  accuracy  regions  for 
two  selected  accuracies;  in  the  latter  case,  the  boundaxy  of  the  stability  region  is  virtually 
indistinguishable  from  the  imaginary  axis.  Each  of  the  figures  carries  a  legend,  specifying 
detailed  stability  characteristics  of  the  scheme.  All  of  the  EuComb'^  schemes  are  A(a)- 
stable,  and  the  legends  specify  the  approximate  values  of  a,  obtained  numerically.  Since  all 
of  the  schemes  EuComb^  are  T-stable  (see  Theorem  4.3),  we  do  not  specify  the  value  of  p 
for  each  of  the  schemes. 

Several  observations  can  be  made  from  Figures  13  —  24. 

1.  There  exist  A-stable  Eulmp^  schemes  of  order  up  to  five  (such  as  EuC'omfel’g).  For 

all  orders  we  have  tested  (up  to  30  or  so),  there  exist  A.(Q:)-stabie  schemes  with  a 
very  close  to  90°.  for  example,  has  order  12  (see  Theorem  4.1),  and 

is  A(Q;)-stable  with  a  >  89.99°.  EuComl^]^  has  order  19  and  is  A(Q;)-stable  with 
a  >  89.996°.  On  the  other  hand,  we  are  not  able  to  predict  reliably  which  of  the 
schemes  will  have  good  stability  properties,  and  which  will  not,  until  such  properties  are 
established  numerically.  EuComh^^^,  for  example,  has  a  >  89.99°,  while  EuCamh^^'^ 
has  89.01°  <  a  <  89.02°.  At  the  other  extreme  is  the  scheme  EuComb^’^,  with  a 
disastrous  a  <  56°  (Figure  21).  As  a  general  rule,  we  have  observed  that  the  schemes 
EuC tend  to  have  poor  stability  properties  whenever  mi+m2  is  even.  This  is, 
of  course,  mostly  a  curiosity,  since  good  schemes  of  various  orders  are  easily  available. 

2.  The  accuracy  regions  for  all  methods  we  tested  are  very  satisfactory,  for  both  real  and 
complex  A.  It  is  easy  to  see  from  Figure  24,  for  example,  that  the  nineteenth  order 
scheme  EuComblo’^g  requires  about  18  nodes  per  wavelength  to  obtain  10  digits. 

6  Adaptive  Implementation 

In  practical  applications  involving  the  numerical  solution  of  ODEs,  issues  such  as  adaptive 
marching,  accuracy  control,  etc.  play  an  important  role.  Since  the  schemes  of  this  paper 
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axe  essentially  of  the  single-step  variety,  naost  issues  arising  in  their  adaptive  implementation 
are  relatively  simple.  It  is  also  worth  keeping  in  mind  that  accuracy  control  is  much  simpler 
when  the  underlying  solver  has  high  convergence  order. 

We  have  implemented  fully  adaptive  versions  of  the  schemes  EuExp,  Eulmp,  EuComb. 
In  this  section,  we  describe  some  of  the  technical  details  of  the  implementation,  while  the 
following  section  describes  some  of  the  numerical  experiments  we  have  performed.  We  will 
be  discussing  these  issues  using  EuExp  as  our  model;  the  other  two  schemes  {Eulmp  and 
EuComb)  encounter  identical  problems,  which  are  handled  in  a  similar  manner. 

6.1  Accuracy  control 

Given  an  initial  value  problem  (1),  (2)  on  the  interval  [a,  6],  an  approximate  solution 
EuExp'^{F,  (f>a),  and  a  positive  e,  we  would  like  to  determine  whether 

I  EuExp^{F, ipa)-ip\<e.  (48) 

Obviously,  this  cannot  be  done  with  complete  reliability;  the  purpose  of  all  existing  tech¬ 
niques  is  to  make  the  determination  with  very  high  probability,  at  an  acceptable  cost.  Fortu¬ 
nately,  the  internal  structure  of  the  method  EuExp^  provides  us  with  a  number  of  conditions 
which  can  be  checked.  Taken  together,  the  conditions  listed  below  have  been  completely  re¬ 
liable  in  our  experience. 

1.  We  verify  that  the  correction  process  has  converged  to  the  precision  e.  In  other  words, 
we  require  that  the  norm  of  the  vector  8  (see  (21)),  obtained  during  the  last  correction, 
be  less  than  e.  This  does  not  guarantee  (48);  it  does  indicate  that  the  correction  scheme 
is  internally  consistent  to  precision  e. 

2.  The  approximate  solution  EuExp'^{F,  (pa)  is  obtained  at  Gaussian  nodes  on  the  interval 
[a,  b].  We  apply  the  operator  W™  to  EuExp'^{F,  (pa),  obtaining  the  m  coefficients  of  its 
Legendre  expansion  (see  (33)).  If  the  discretization  of  EuExp^{F,<pa)  is  sufficiently 
fine,  the  last  several  coefficients  of  the  Legendre  expansion  must  be  small.  In  our 
implementation,  we  demand  that  the  last  two  coefficients  be  smaller  than  e. 

3.  Yet  another  test  we  perform  attempts  to  verify  simultaneously  that  both  the  correction 
process  and  the  discretization  have  converged  to  precision  e.  Once  EuExp'^{F,(pa) 
has  been  evaluated,  we  obtain  the  value  of  the  approximate  solution  at  the  point  b 
via  interpolation  (see  Observation  4.1).  We  apply  the  interpolation  process  to  both 
EuExp^{F,(Pa)  and  EuExp^^{F,p>a)i  demand  that  the  difference  be  less  that  e. 

4.  In  extreme  cases,  the  solution  of  the  ODE  can  be  so  underresolved  as  to  become 
unstable;  the  usual  result  is  exponential  overflow.  To  guard  against  this  condition. 
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we  check  the  size  of  the  solution  at  every  step  of  the  Euler  process;  if  the  solution  is 
sufficiently  large  (we  have  arbitrarily  set  the  threshold  to  10^®)  the  point  of  view  is 
taken  that  the  problem  is  underresolved. 

6.2  Step- length  control 

Our  approach  here  is  completely  standard.  We  start  with  a  more  or  less  arbitrary  step-size, 
and  attempt  to  apply  the  scheme  EuExp-^.  If  the  resulting  precision  is  insufficient  (according 
to  the  criteria  (a)  -  (d)  above),  the  step-length  is  halved,  and  the  process  repeated.  If  the 
precision  is  satisfactory,  the  steplength  is  unchanged.  If  the  precision  is  satisfactory  two 
steps  in  a  row,  the  step-length  is  doubled. 

6.3  Linearly  implicit  implementation 

In  order  to  reduce  the  number  of  function  evaluations,  a  common  practice  in  most  extrapo¬ 
lation  codes  is  to  use  a  “linearly  implicit”  formulation  of  the  marching  scheme  [10].  For  the 
ordinary  differential  equation 

=  Fit,<p{t))  te[a,b]  (49) 

with  an  approximate  solution  (po{t),  we  let  (p  =  ifo  + 6  and  write  (49)  in  the  form 

foit)  +  ^'{t)  =  +  ^{i)) 

or 

S{t)=  f  F{T,(po{T)  +  S{T))dT -(po{t).  (50) 

J  a 

But  for  small  6,  we  have 

+  «(i))  =  +  J„(()<(t)  +  0(||i||2), 

where  denotes  the  Jacobian  of  F  with  respect  to  its  second  argument  at  the  point 

(t,  (/:>o(t)).  Substituting  this  expression  into  (50),  we  get 

6{t)=  f  J^^{t)S{T))dT+  I  F{T,<po(T))dT-(po{t).  (51) 

We  can  then  use  the  backward  Euler  method  to  drive  a  deferred  correction  process  applied 
to  (51).  Our  linearly  implict  code  performs  up  to  six  steps  of  deferred  correction  on  this 
equation,  after  which  the  approximate  solution  (po  is  updated  according  to 

Voit)  :=  (po{t)  +  S{t). 
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Table  1:  Performance  of  low,  moderate,  and  high  order  spectral  deferred  correction  methods 
for  non-stiff  problems.  The  first  column  indicates  the  requested  precision  and  the  remaining 
columns  list  the  number  of  function  calls  required  by  the  corresponding  scheme. 


Precision 

EuExp\ 

EuExpl 

EuExpIq 

10"^ 

70 

44 

93 

10"® 

287 

176 

155 

10-12 

— 

2574 

310 

7  Numerical  Experiments 

We  illustrate  the  performance  of  spectral  deferred  correction  methods  with  two  examples. 
The  first  is  the  system  of  three  ordinary  differential  equations  satisfied  by  the  Jacobian 
elliptic  functions  sn,  cn,  dn: 


sn\t)  —  cn{t)  •  dn{t) 
cn'{t)  =  —sn{t)  •  dn{t) 
dn\t)  =  —n  •  sn{t)  •  cn{t) 

with  fi  =  0.5  on  the  interval  [0,1]  with  initial  data  sn(0)  =  0,  cn(0)  =  1,  dn(0)  =  1.  This 
is  a  common  model  for  nonstiff  problems.  As  can  be  seen  from  the  results  in  Table  1,  the 
order  of  accuracy  of  the  method  should  increase  with  the  desired  precision  to  obtain  optimal 
performance. 

Our  second  example  is  the  Van  der  Pol  oscillator,  a  well-known  stiff  system  of  two  equa¬ 
tions: 


y[{i)  =  y2{t) 

y'iit)  =  {l-yl{t)y2{t)-yi{t))/e 

with  t/i(0)  =  2,  2/2(0)  =  0.  We  choose  e  =  10~®  and  solve  on  the  interval  [0,1].  For  the 
sake  of  comparison,  we  tested  our  codes  against  one  of  the  best-performing  codes  for  this 
problem  -  the  high  order  extrapolation  code  SEULIM,  due  to  Deuflhard,  Nowak  and  Poehle 
[10].  Our  results  are  collected  in  Tables  2-4. 

A  number  of  observations  can  be  made  from  these  tables. 

1.  The  code  SEULIM  requires  noticeable  fewer  function  evaluations  than  the  Eulmp 
schemes  at  any  requested  precision.  SEULIM  also  requires  fewer  function  evaluations 
than  the  linearly  implicit  deferred  correction  code. 
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Table  2:  Performance  of  the  extrapolation  code  SEULIM  on  the  Van  der  Pol  oscillator 
problem.  The  first  column  indicates  the  requested  precision,  the  second  column  indicates 
the  number  of  function  evaluations,  the  third  column  lists  the  computed  solution  component 
2/1(1)  and  the  fourth  column  lists  the  computed  solution  component  2/2(1)* 


Precision 

F.  calls 

J/i(l) 

^2(1) 

10-^ 

166 

-1.79765618619 

0.000017424290 

10"=^ 

250 

-2.00019319529 

0.000060543984 

10-3 

751 

1.98059515302 

-0.033130626096 

lO-'* 

1197 

1.70727095158 

-0.885933089211 

10-5 

1872 

1.70633329346 

-0.892323342916 

10-3 

2622 

1.70618445916 

-0.892674570523 

10-7 

3870 

1.70616982299 

-0.892804713111 

10-3 

5476 

1.70616796672 

-0.892809443466 

10-® 

6531 

1.70616775897 

-0.892809533238 

0 

1 

0 

10490 

1.70616773572 

-0.892809663387 

10-^^ 

21547 

1.70616773312 

-0.892809699177 

10-^2 

86899 

1.70616773227 

-0.892809700614 

Table  3:  Performance  of  the  spectral  deferred  correction  code  Eulmp  on  the  Van  der  Pol  os¬ 
cillator  problem.  The  first  four  columns  correspond  to  those  in  Table  2.  The  last  two  columns 
show  the  number  of  points  n  used  on  each  subinterval  (the  maximal  order  of  accuracy)  and 
the  maximal  number  of  corrections  ncorr  actually  used  by  the  code. 


Precision 

F.  calls 

J/i(l) 

J/2(l) 

n 

ncorr 

10-^ 

3462 

1.70645342840 

-0.892494581060 

6 

5 

10-2 

5340 

1.70642671028 

-0.892530184336 

8 

4 

10-3 

6754 

1.70615323588 

-0.892825220721 

8 

4 

10-^ 

11920 

1.70616745993 

-0.892809992650 

16 

8 

10-5 

17880 

1.70616771499 

-0.892809719348 

22 

12 

10-3 

20576 

1.70616773089 

-0.892809702550 

22 

12 

10-" 

21852 

1.70616773187 

-0.892809701498 

22 

12 

10-3 

23366 

1.70616773221 

-0.892809701142 

22 

12 

10-^ 

29798 

1.70616773217 

-0.892809701012 

22 

12 

10-10 

56562 

1.70616773217 

-0.892809701047 

22 

12 

10-“ 

128362 

1.70616773217 

-0.892809701031 

22 

12 
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Table  4:  Performance  of  the  linearly  implicit  spectral  deferred  correction  code  on  the  Van 
der  Pol  oscillator  problem.  The  columns  correspond  to  those  in  Table  3. 


Precision 

F.  calls 

yi(l) 

2/2(1) 

n 

ncorr 

10-^ 

1976 

1.70756400695 

-0.891306260408 

6 

5 

10-2 

2587 

1.70629621863 

-0.892672011388 

8 

4 

10-3 

3112 

1.70618777664 

-0.892788115757 

8 

4 

10-“ 

4203 

1.70616787726 

-0.892809527001 

16 

8 

10"® 

4839 

1.70616773785 

-0.892809694425 

22 

12 

10-® 

5644 

1.70616773478 

-0.892809697543 

22 

12 

10-" 

5887 

1.70616773216 

-0.892809699781 

22 

12 

10-® 

7817 

1.70616773207 

-0.892809701150 

22 

12 

0> 

1 

o 

13885 

1.70616773217 

-0.892809700662 

22 

12 

10-10 

28857 

1.70616773217 

-0.892809700900 

22 

12 

10-“ 

54528 

1.70616773217 

-0.892809701003 

30 

20 

2.  If  actual  accuracies  are  compared,  rather  than  requested  precision,  a  slightly  different 
picture  begins  to  emerge.  The  deferred  correction  scheme  achieves  about  eight  digits  of 
accuracy  at  a  requested  tolerance  of  10“®,  using  4,839  function  calls.  SEULIM,  on  the 
other  hand,  acheives  eight  digits  of  accuracy  at  a  requested  tolerance  of  10“^°,  using 
10,490  function  calls.  This  is  not  intended  to  disparage  SEULIM’s  performance.  Our 
implementation  simply  has  very  strict  (perhaps  excessive)  error  control.  If  CPU  times 
were  compared,  SEULIM  would  be  the  clear  winner. 

8  Conclusions 

We  believe  that  deferred  correction  methods  based  on  an  integral  equation  formulation  of 
the  ordinary  dijfferential  equation  are  promising  candidates  for  further  investigation.  They 
have  excellent  stability  properties,  are  easy  to  implement,  and  require  only  a  good  low- 
order  solver  to  drive  the  process.  Our  preliminary  experiments,  using  a  primitive  adaptive 
implementation,  compare  favorably  with  a  state-of-the-art  extrapolation  code  at  moderate 
to  high  precision. 
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Stability  and  accuracy  regions  for  Detail  of  stability  and  accuracy 

Eulmp^',  fj,  ^  —.3913,  a  =  90®  regions  for  Eulmp\ 


Figure  7:  Figure  8: 

Stability  and  accuracy  regions  for  Detail  of  stability  and  accuracy 

EuImpQ]  (X  ^  —.3101,  a  ^  89.979®  regions  for  Eulmp\ 


Figure  13: 

Figure  14: 

Stability  and  accuracy  regions  for 

Detail  of  stability  and  accuracy 

EuComb^’^’,  a  =  90° 

regions  for  EuCombQ^ 
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Figure  15: 

Stability  and  accuracy  regions  for 
EuComb\l’^2i  ^  ^  89.9914° 


Figure  16: 

Detail  of  stability  and  accuracy 
regions  for  EuCombi^’^^2 


Figure  19:  Figure  20: 

Stability  and  accuracy  regions  for  Detail  of  stability  and  accuracy 

EuComb\^’^Q]  a  «  89.014°  regions  for  EuComb\^’^g 


Figure  23: 

Stability  and  accuracy  regions  for 
EuC omblQ’^l]  a  w  89.9969® 


Figure  24: 

Detail  of  stability  and  accuracy 
regions  for  EuComb^'^g 


